/*

  Cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - No memory coalescing in the Row reads
  - Coalescing memory writing through transposing in shared memory
  - All threads read in all data

 */

CUBE_KERNEL(CoalescingTriangle, cuFloatComplex *array_d, cuFloatComplex *product_d)
{
  CUBE_START;

  //set global thread ID
  if ((blockIdx.x)*TILE_WIDTH > (blockIdx.y)*TILE_HEIGHT) return;
  int Row = blockIdx.y*TILE_HEIGHT + threadIdx.y ;
  int Col = blockIdx.x*TILE_WIDTH + threadIdx.x;

  //declare shared memory for output
  __shared__ float4 output[TILE_WIDTH*TILE_HEIGHT*NPOL];
  //get local thread ID
  int localRow = threadIdx.y;
  int localCol = threadIdx.x;

  //instantiate sum variables
  cuFloatComplex sumXX = make_cuFloatComplex(0.0, 0.0);
  cuFloatComplex sumXY = make_cuFloatComplex(0.0, 0.0);
  cuFloatComplex sumYX = make_cuFloatComplex(0.0, 0.0);
  cuFloatComplex sumYY = make_cuFloatComplex(0.0, 0.0);
  
  for(int t=0; t<NTIME; t++){
    float4 blurbRow = ((float4*)array_d)[t*NSTATION + Row];
    float4 blurbCol = ((float4*)array_d)[t*NSTATION + Col];
    cuFloatComplex rowX = make_cuFloatComplex(blurbRow.x, blurbRow.y);
    cuFloatComplex rowY = make_cuFloatComplex(blurbRow.z, blurbRow.w);
    cuFloatComplex colX = make_cuFloatComplex(blurbCol.x, blurbCol.y);
    cuFloatComplex colY = make_cuFloatComplex(blurbCol.z, blurbCol.w);
    sumXX = cuCaddf(cuCmulf(rowX, cuConjf(colX)), sumXX);
    sumXY = cuCaddf(cuCmulf(rowX, cuConjf(colY)), sumXY);
    sumYX = cuCaddf(cuCmulf(rowY, cuConjf(colX)), sumYX);
    sumYY = cuCaddf(cuCmulf(rowY, cuConjf(colY)), sumYY);
  } 


  output[(localRow*TILE_WIDTH + localCol)*NPOL + 0] = make_float4(sumXX.x, sumXX.y, sumXY.x, sumXY.y);
  output[(localRow*TILE_WIDTH + localCol)*NPOL + 1] = make_float4(sumYX.x, sumYX.y, sumYY.x, sumYY.y);

  __syncthreads();

  int rowIndex = Row*(Row+1)/2;
  for (int pol=0; pol<NPOL; pol++) {
    // new mapping of threads to columns
    localCol = (pol*TILE_WIDTH + threadIdx.x)/2;

    // need to check if the column is greater than the row
    if ((blockIdx.x)*TILE_WIDTH + localCol > Row) return;

    ((float4*)product_d)[(rowIndex*NPOL + (pol + blockIdx.x*NPOL)*TILE_WIDTH) + threadIdx.x] =
      output[(localRow*NPOL + pol)*TILE_WIDTH + threadIdx.x];
  }

  CUBE_ADD_FLOPS(NTIME*NPOL*NPOL*8);
  CUBE_ADD_BYTES(NTIME*NPOL*8 + NPOL*NPOL*8);

  CUBE_END;
}
